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Chapter 1 

Introduction 


Unlike the uniform density spherical shell appoximations of Newton, the con- 
sequence of spaceflight in the real universe is that gravitational fields are sen- 
sitive to the nonsphericity of their generating central bodies. The gravitational 
potential of a nonsplrerical central body is typically resolved using spherical 
harmonic approximations. However, attempting to directly calculate the spher- 
ical harmonic approximations results in at least two singularities which must be 
removed in order to generalize the method and solve for any possible orbit, in- 
cluding polar orbits. Three unique algorithms have been developed to eliminate 
these singularities by Samuel Pines [1], Bill Lear [2], and Robert Gottlieb [3]. 

This paper documents the methodical normalization of two 1 of the three 
known formulations for singularity-free gravitational acceleration (namely, the 
Lear [2] and Gottlieb [3] algorithms) and formulates a general method for defin- 
ing normalization parameters used to generate normalized Legendre Polyno- 
mials and ALFs for any algorithm. A treatment of the conventional formula- 
tion of the gravitational potential and acceleration is also provided, in addition 
to a brief overview of the philosophical differences between the three known 
singularity- free algorithms. 


1 The Pines algorithm (Section 4.1) lias been previously normalized and thoroughly in- 
vestigated by Lundberg and Schutz [4] and subsequently implemented by DeMars [5]. See 
Section 4.1.2. 
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CHAPTER 1. INTRODUCTION 
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Chapter 2 

Gravitational Potential and 
Acceleration 

2.1 Coordinates 

Gravitational potential V is typically resolved using an equatorial central-body- 
fixed Cartesian coordinate system [ xi, yb Zb \ , where 

• Xb passes through the center of mass, equator, and prime meridian of the 
central body 

• Zb is the north polar axis of the central body 

• yb completes the right-handed triad, where yb is positive in the eastern 
hemisphere 

Spherical coordinates [ r 9 <f> ] T can then be defined in terms of the 

central-body-fixed coordinates 

• r = \J x\ + yl + z% = magnitude of position vector 

• 9 = arctan2(j/h, Xb) = (positive east) longitude or right ascension 

• <f> = arcsin^) = (positive north) latitude or declination 

If Xb = yb = 0, then 9 may be set to any value. For convenience, it is 
typically set to 9 = 0 in this case. Similarly, if r = 0, then 0 may be set to any 
value but is typically set to (j> = 0. 

The relationship between central-body-fixed and spherical coordinates is 
shown in Figure 2.1. The spherical coordinate angles themsevles are rarely used 
directly since they typically appear as arguments in trigonometric functions. 
See Section 2.1.2 for more information. 
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4 CHAPTER 2. GRAVITATIONAL POTENTIAL AND ACCELERATION 


Zb 



Xb 


Figure 2.1: Body-fixed coordinates and spherical coordinates illustrated 

2.1.1 The Gravitational Acceleration 

The gradient of the gravitational potential V is the force of gravity per unit 
mass, which is the acceleration. 

The gravitational acceleration in the central-body-fixed Cartesian coordinate 
system [ x b yb z b ] is 

r T w T/ r dv dV dV 1 T /o i \ 

a — [ a Xb a Vb a Zb J — V b V — ^ dxb dyb dzb J (2.1) 

where V& is the gradient with respect to the central-body- fixed coordinates. 

Taking the gradient with respect to a central-body-fixed (i.e. accelerating) 
coordinate system means that we must include the Coriolis, centripetal, tangen- 
tial, and relative accleration terms [6, pg. 54-55]. To eliminate the complication 
of including these terms, an alternative inertial (non-accelerating) coordinate 
system centered in the central body is used to compute the gradient. 

Consider a displacement Ad of the position vector in spherical coordinates 
in the direction of increasing 9 (away from the x b axis). While displacement of 
9 is actually the arc of a circle, a limit is approached as A 9 — > 0 when taking the 
gradient in spherical coordinates and thus the displacement arc becomes arbi- 
trarily close to its subtending chord. The instantaneous positive displacement of 
9 for a vector in spherical coordinates is thus perpendicular to the radial vector 
in the direction of increasing 9. This logic can be similarly applied for (j). The 
directions of increasing 9 and (f> form the basis of the new coordinate system [7, 
pg. 143]. 

The new x axis is parallel to the position vector and is defined as the r 0 
axis. The new y axis is located on the central-body-fixed Xb‘y b plane at a right- 
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2.1. COORDINATES 


5 



Figure 2.2: Orthogonal spherical coordinates [ r a 8 a <j) 0 ] T and central- 
body-fixed coordinates [ Xb yb Zb ] T shown with the r 0 (f> 0 plane 

angle to the r a axis in the direction of increasing 8 and is defined as the 9 0 
axis. The new 2 axis completes the right-handed triad and is defined as the 
4> 0 axis. The (j) 0 axis forms the angle </> with the zi, axis in the direction of 
increasing (f) (away from the Zb axis). The new coordinates are referred to as 
the orthogonal spherical coordinates [ r 0 8 0 (f> 0 ] , not to be confused with 

the conventional spherical coordinates [r 8 (f> ] T . The relationship between 
central-body-fixed coordinates and orthogonal spherical coordinates is shown in 
Figure 2.2. 

The transformation matrix from the orthogonal spherical coordinates to the 
central-body-fixed coordinates is [2] 

cos cf) cos 8 — sin 9 — sin (j> cos 8 r 0 

cos(j)S\n8 cos 9 — sinc(>sin0 9 0 (2-2) 

sin 4> 0 cos (j) 4>o 

All of the trigonometric functions in the transformation can be precomputed 
using the trigonometric relationships outlined in Section 2.1.2. 

The gravitational acceleration in the orthogonal spherical coordinate sys- 
tem is [6, pg. 54] 

a — [ 0-r o a,e o a r f >0 J — v V — | g r rcosj> d6 r d<t> J (2->) 

The result of Equation 2.3 can then be transformed back to central-body-fixed 
coordinates via Equation 2.2. 
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6 CHAPTER 2. GRAVITATIONAL POTENTIAL AND ACCELERATION 


2.1.2 Trigonometric Relationships 

Instead of solving for the spherical coordinate angles 4> and 9 , it is more efficient 
to precompute the values of trigonometric functions of these angles given by the 
definitions of sine and cosine. 


sin 9 = 

Vb 

(2.4) 

V x b + Vb 

cos 9 = 

x b 

(2.5) 

V x b + y'b 

sin (f> = 

Zb 

r 

(2.6) 

COS (f) = 

V x b + y'b 

(2.7) 


A logical check should be included prior to computing these values to identify 
the cases Xb = yb = 0 or r = 0. If Xb = yb = 0, then the equations sin0 = 0 
and cos 9 = 1 are typically substituted. Similarly, if r = 0, then the equations 
sin^ = 0 and cos <f> = 1 are typically substituted. 


2.2 Constants and Coefficients 

S n ,m and C n , m are known in the literature as the spherical harmonic mass coeffi- 
cients of the central body, which we will refer to simply as the mass coefficients. 
Certain subsets of the mass coeffients are given special names. 

• Cn,o = zonal coefficients 

• Sn,mC n)n = sectorial coefficients 

• Sn,m,C n)rn = tesseral coefficients 

Tables of mass coefficients are usually provided with the gravitational pa- 
rameter /i = CM (where G is the Newtonian gravitational constant and M is 
the mass of the central body) and the scaling radius a eq for which the coefficients 
are calibrated. 


2.3 Legendre polynomials and ALFs 

The various functions denoted P in this paper are the Legendre polynomials and 
the Associated Legendre Functions (ALFs). Legendre polynomials of the first 
kind are denoted P n with argument sin </>. ALFs of the first kind are similarly 
denoted P n ,m also with argument smtj). Legendre polynomials and ALFs have 
subscipts n and m which are called the degree and order of the polynomial, 
respectively. ALFs with m > n are defined as zero. The Legendre polynomials 
are equivalent to the ALF of the same degree but with m = 0. 

In the remainder of this document, the argument of sin 4> will be assumed 
and omitted for brevity for all Legendre polynomials and ALFs. 
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2.4 The Gravitational Potential Function 

The gravitational potential V satisfies Laplace’s equation 

d 2 V d 2 V d 2 V n 

9 T o T 9 — 0 

dx b 2 dy 2 dz b 2 


(2.8) 


The potential function V can be written in spherical coordinates as an or- 
thogonal expansion using spherical harmonics [6, Pg. 52]. 




oo n 


1 + EE(v) > " SIP TTlO Cri^m COS TTlO^j 


n— 1 m — 0 


(2.9) 


Equation 2.9 can also be written with the zonal terms separated out 


v = a 


£(v) 


^ p n c, 


nW, 0 


n— 1 
oo n 


+ ££(v)”^ ,m ( Sn,m sin md + C n ^ m cos mO) 


n = 1 m = 1 


( 2 . 10 ) 


2.4.1 Square Gravity Models 

The S ntm , C nrn . and P n>m values can be used to form lower triangular matricies, 
where the terms above the diagonal are all zero, n is the row, and m is the 
column. 

When the origin of the central body-fixed coordinate system is located at the 
center of mass, the C qo, Si,i, and Cyi coefficients are all zero. This convention 
will be adopted throughout this paper. Therefore, in computing V, the sums 
may be made starting with n = 2. Models where maximum degree and order 
are equal are referred to as square models. Let rid be the desired maximum 
degree and order of the gravity model. Then 


V = 




1 + E ( ( Sn ’ m sin m 9 + Cn ' m cos m 

n— 2 m = 0 

With the zonal terms separated out, Equation 2.11 becomes 

i 1 + E(v)” p ” c ”'“ 

n — 2 

££(v)" p " ,m ( S n ,m sin m6 + C n , m cos md) 


( 2 . 11 ) 




n — 2 m=l 


( 2 . 12 ) 
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2.4.2 Non-Square Gravity Models 

Whereas Equations 2.11 and 2.12 assume square models, it is also possible to 
approximate the gravitational potential using non-square models. This can be 
accomplished by changing the bounds of the sums in Equations 2.11 and 2.12 
to nd and m,j, the desired degree and desired order, respectively. The outcome 
of a non-square model can then be written as 


V= £ 


rid ri 

i + E E ,771 (^ 77,777 sin TflO H - Cri,m COS TYlO ) 


77=2 777=0 
777 < 777 ( i 


resembling Equation 2.11, or separating out zonal terms, 


v = ^ 


nd 


E(^) 


77=2 

rid n 


+ E (~) Pn,m(Sn,mSmmd + C n , m cosmO) 


77 = 2 777=1 
777 < 777 


(2.13) 


(2.14) 


resembling Equation 2.12. 

Because most algorithms for generating the Legendre polynomials and ALFs 
are optimized to generate values for square models, it is best to pass only the 
desired degree of the non-square gravity model to these function-generating 
subroutines. In calculating a non-square potential, the excess ALFs beyond the 
desired order are then simply unused. 

The fact that Equations 2.11 and 2.12 (and thus Equations 2.13 and 2.14) 
are orthogonal expansions of V means that any lower order expansion is merely 
a truncated higher order expansion. No refit of the coefficients is necessary [2]. 
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Chapter 3 


Normalization 


3.1 The Normalization Factor 

Using a normalization factor allows mass coefficients to be electronically rep- 
resented as tractable values well within the valid range for IEEE 754 double- 
precision floating point variables, even when degree n and order m are relatively 
large. The mass coefficients of a central body are “normalized” when they are 
divided by this normalization factor iV„ jTO , typically defined as [8, pg. 544] 


N = 


'{n - m)\(2n + 1)(2 - <5 0 , m ) 
(n + to)! 


N n = N, 


n, 0 


(3-1) 


where So,m is the Kronecker delta function that returns one if to = 0 and zero 
otherwise. Table 3.1 lists the normalization factors through degree and order of 
four. 

Historically, mass coefficients were “unnormalized” by the end user of the 
model by multiplying the mass coefficients with their corresponding normaliza- 
tion factor. This was done because conventional gravitational potential algo- 
rithms required unnornralized coefficients in their formulations. A fundamental 



TO — > 

n 

0 

1 

2 

3 

4 

0 

1 

0 

0 

0 

0 

1 

73 

Vs 

0 

0 

0 

2 

75 

/? 
V 3 

1 / 3 " 
2 V 3 

0 

0 

3 

77 

[7 

V 6 

1 r 

2 y 15 

1 

6V 10 

0 

4 

3 

3 \[T. 

1 fi 

1 r~^~ 

1 r~^~ 



10 

2 y 5 

2 V 70 

8 V 35 


Table 3.1: Normalization factors through 4x4 
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CHAPTER 3. NORMALIZATION 


problem encountered by attempting to unnormalize mass coefficients with very 
high degrees and orders is that the normalization factors become prone to over- 
flow in modern computers. This is due to the factorial in the denominator of 
Equation 3.1, (n + m)\. The user is restricted to n ( j + m.d < 171 using IEEE 
754 double-precision floating point variables to calculate normalization factors 
because of overflow beyond this boundary (see Section 3.2). 

It was later established that the additional work of unnormalizing mass coef- 
ficients becomes unnecessary by “normalizing” ALFs in gravitational potential 
formulations. Legendre polynomials and ALFs are said to be “normalized” when 
multiplied by N n ^ m . This normalization scheme is ideal because the product of 
an ALF and its corresponding coefficient is equal to the product of the normal- 
ized ALF and its corresponding normalized coefficient. In this paper, an overbar 
will be used to indicate normalized quantities, i.e. C n rn and P n , m - By the afore- 
mentioned normalization conventions, Equation 3.2 shows that the product of 
the ALF and corresponding coefficient holds true for normalized quantities. 

Pn,mCn,m — 

This eliminates the restrictions imposed by requiring unnormalized mass coef- 
ficients for gravitational potential formulations. The relationship from Equa- 
tion 3.2 can be seen with S n>m and S n>m as well. 



3.2 Recursive Mass Coefficient Normalization 

Gravity models evaluating spherical harmonic associated Legendre functions 
at high degree n and order m require normalized Legendre coefficients to ac- 
curately compute terms ranging over hundreds of orders of magnitude. The 
standard normalization scheme is to divide the unnornralized coefficients by the 
Kaula normalization factor N n>m , defined in Equation 3.1. If N n , m is calcu- 
lated directly, computation difficulties arise when n + m > 170 because 171! 
triggers an overflow condition in IEEE double precision (64-bit) real numbers. 
The overflow limit is ±1.797693134862316e+308 1 . 

Suppose N n>m values are directly calculated and grouped as elements in a 
lower-triangular matrix, as is customary, by incrementing m from 0 until in = n 
before n is incremented to begin another matrix row. In performing this task, 
the first ( n + m)\ overflow will be encountered for N 8885 . Table 3.2 provides 
N n rn values neighboring the N 8 6,85 element. In Table 3.2, N n>rn values which can 
be computed are nowhere near an overflow condition because they are quotients 
with large denominators. This suggests a recursive computation will succeed in 
generating accurate N n ^ m values far beyond element N 8 6 ^ 5 . The recursion will 
only fail when the N n<m quotient overflows. 

To document the recursion, suppose the value of N n i trn / is given as x. As- 

1 This number was obtained from the realmax( ‘double ’ ) command in MATLAB. 
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Degree n 

84 

Order m 
85 

86 

85 

1.117258027e+151 

1.456726244e+152 

N/A 

86 

1.024089587e+152 

Overflow 

Overflow 

87 

Overflow 

Overflow 

Overflow 


Table 3.2: Normalization factors neighboring A^ 6,85 





Order 



Degree 


ml — 1 

m! 


m' + 1 

n' — 1 


j (n'— m')(2n / +l) 

X \J (n' -\-m' )(2n' — 1) 


n' 

X M 

/ ^ <$o ,m' 

X 

X \J 

/ (n / +m / +l)(n / — m') 

(n / +m / )(n / — m' + l) 


n! + 1 


/(n' + l+m')(2n' + l) 
V (n' + l-m')(2n'+3) 



Table 3.3: Recursions for generating normalization factors 


suming n! — m! > 0 for “upward” or “rightward” recursion 2 , Table 3.3 supplies 
recursive formulae for adjacent elements in terms of x and its associated de- 
gree n' and order in' . Each formula has been verified using the three finite 
N n m elements appearing in Table 3.2. For example, invoking the “downward” 
recursion, 

-^86,84 = -^85,84^/- ^ ~ ~ 9. 16609737241 _/V 8 5, 84 

Most gravity models are formally published with normalized coefficients. 
Those wishing to use normalized coefficients at n + m > 170 in an unnormalized 
model will find the Table 3.3 recursions useful in their work. At n+m < 170, the 
Table 3.3 recursions offer computational efficiencies over Equation 3.1 evalua- 
tions, but these will be of little consequence if coefficients are to be unnormalized 
in a single pass with the results stored for all subsequent use. 


3.3 Normalization Ratios 

Consider a recursion formula for Legendre Polynomials [9, pg. 114, sec. 3] 

P n = ^[(2n - 1) sin(/>P„_i - (n - 1)P„_ 2 ] (3.3) 

The corresponding normalized Legendre polynomial P n is found by multiplying 
Equation 3.3 by the normalization factor N n . 

P n = N n P n = -[(2 n - 1) sirup N n P n _i - (n - l)V„P n _ 2 ] (3.4) 

n 

2 These forbidden recursions would otherwise step outside lower-triangular matrix limits if 
x corresponded to a diagonal element with n' — m! = 0. 
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CHAPTER 3. NORMALIZATION 


Because Equation 3.4 returns normalized polynomials and must recur over 
its own inputs, the equation needs to be written as a function of normalized 
polynomials by replacing the conventional polynomials with their normalized 
equivalents. By substituting 


Equation 3.4 becomes 


p _ 

n N n 


(3.5) 


Pn = 


n 


N - N 

(2 n- 1) sin (f> n P n _ ± - ( n - 1) n P ra _ 2 

Wn-1 ^n - 2 


(3.6) 


Equation 3.6 now contains ratios of normalization factors, namely and 

Nr, 


3.4 The Recursion Normalization Parameter 


The ratios of the normalization factors, which are referred to as the normal- 
ization parameters and denoted by A, can be pre-computed since many ratios 
will be needed more than once while normalizing the algorithms. The ratios of 
Equation 3.6 are written as A„_i and A„_ 2 , where the subscripts of the parame- 
ters refer to the subscript of the corresponding polynomial which they normalize 
in Equation 3.6. 

When n = 2, A„_i = X\ . However, when n — 3, A „_2 ^ X\- For this reason, 
the parameters are written as functions of the current values of n and m , such 
as A„_i(n) and \ n - 2 , m (n, m), and ignore the actual value of the subscripts of 
the parameters. The subscripts are merely notation used to identify a specific 
parameter for normalizing the ALF with the same subscripts. 

This notation is used in Equation 3.6 to obtain 

P n = -[(2 n- l)sin^A n _i(n)P n _i - (n - 1)A„_ 2 (n)P„_ 2 ] (3.7) 

n 

An equation for the parameter A„_i can be found by substituting the defi- 
nition of the normalization factors (Equation 3.1) in the ratio and simplifying. 


N n \/2n + 1 l2n+l 

~ N n _ i “ ^/2 (n - 1) + 1 _ V 2 n - 1 


(3.8) 


The Lear algorithm (Section 4.2) requires five normalization parameters to 
recursively compute normalized Legendre polynomials and ALFs. Each of these 
parameters can similarly be derived by simplifying the definition of the normal- 
ization factor for each of the corresponding ratios. 
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Chapter 4 

The Three Singularity-Free 
Algorithms 


By inspection, we can identify two potential singularities in Equation 2.3. They 
can occur when cos <f> = 0 (such as in a polar orbit) or when taking the partial 
derivative of ALFs with m = 1. Three unique algorithms have been developed 
to eliminate these singularities by Samuel Pines [1], Bill Lear [2], and Robert 
Gottlieb [3]. 

4.1 Pines Algorithm 

4.1.1 Basis of the Pines Approach 

Pines approached the problem by first transforming a position vector to central- 
body-fixed coordinates. By reallocating factors in each term of the potential, 
he defined a series of special functions in terms of the unit position vector 
components which can then be solved recursively without singularities. A set 
of polynomials he referred to as the derived ALFs were created by modifying 
the conventional definition of ALFs. The derived ALFs have similar recursive 
behaviors to their conventional counterparts but have no discontinuities in their 
partial derivatives. 

4.1.2 Pines Algorithm Implementations 

The Pines acceleration algorithm has previously been normalized by Lundberg 
and Schutz [4]. Within the normalized algorithm, Pines’ derived ALFs can be 
recursively generated in a number of ways. To evaluate the various approaches, 
Lundberg and Schutz developed seven different recursions algorithms and then 
performed numerical analyses of the stability of normalized and unnormalized 
versions of each recursion algorithm. Lundberg and Schutz concluded that a 
simple row-wise or column-wise recursion provides the most stability of all the 

13 
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,„_i(n) = 

N n _ [ 

c 

,„_ 2 (n) = 

N n _ [ 

1 

to 

1 

< 

,n-l (n) = 

N n , n 

-**71—1,71—1 


^ (n, m) = 

N 

1 y n,m 

N i 

1 y n— l,m 


= 

N 

1 y n,m 

N n,— 2,ra 


2n + 1 
2n — 1 
2n -f- 1 
2n — 3 


1 /2n + l 


2n — 1 V 2n 


(n — m)(2n + 1) 


(n + m)(2n — 1) 

I (n — m)(n — m — 1)(2 n + 1) 


(n + m)(n + m— 1)(2 n — 3) 

Table 4.1: Normalization parameters A for a normalized Lear algorithm 


(4.1) 

(4.2) 

(4.3) 

(4.4) 

(4.5) 


algorithms, normalized or unnormalized. In this analysis, the normalized and 
unnormalized versions of the column-wise recursion by Lundberg and Schutz [4, 
Recursion I] are utilized for implementations of Pines. The implementations 
used here are based on a normalized implementation provided by DeMars [5]. 


4.2 Lear Algorithm 

4.2.1 Basis of the Lear Approach 

Lear transformed the position vector to the orthogonal spherical coordinate 
system. This results in several sec (f) factors that emerge in the equations for 
the acceleration. Lear found that stable recursions could be developed by as- 
similating the sec (j> factors in the ALF recursion equations. Lear utilized the 
conventional (unmodified) ALFs and thus used traditional recursion equations 
for ALFs, across which he simply distributed the sec <j> factors. Terms in the 
potential which include an ALF but no sec (f> could easily eliminate the secant 
(which has been combined in the value of the ALF) by multiplying the term by 
cos </>, a value with no discontinuity. 

4.2.2 Normalized Lear Algorithm 

Each of the recursion equations of the Lear algorithm is now normalized using 
the normalization parameters of Section 3.4. The parameters required in the 
Lear algorithm are defined in Table 4.1. Each of the original equations can be 
found in Ref. [2]. Only the normalized equations are presented here and should 
replace their analogues in the original algorithm. The added normalization 
parameters are denoted with an underbrace or overbrace to bring attention to 
the changes from the equations in the original document. 
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In each of the following recursions, let rid be the desired degree and order of 
the gravity model. 

Recursion for Zonal Legendre Polynomials P n 

This recursion utilizes Parameters 4.1 and 4.2. For n = 2 through nd 

P n = 1) sin <j> A„_i(n)P„_i - (n - 1) A n _ 2 (n) P„_ 2 ] (4.6) 

4.1 4.2 

where Po = 1 and Pi = N\P\ = y/3 sin<^. 

Recursion for Zonal Legendre Polynomial Derivatives P' n 

This recursion utilizes Parameter 4.1. For n = 2 through nd 

P'n = A„_i(n)[sin^P'„_i + nP„_i] (4.7) 

4.1 

where P'i = \/3. 

Recursion for Sectorial ALFs (sec 4>P n , n ) 

This recursion utilizes Parameter 4.3. For n = 2 through nd 

(sec (j)Pn, n ) = (2 n- 1) cos (j) A„_ li „_ 1 (n)(sec^P„_ 1>n _ 1 ) (4.8) 

4.3 

where (sec</>Pi,i) = V^- 

Recursion for Tesseral ALFs (sec0P„ jm ) 

This recursion utilizes Parameters 4.4 and 4.5. For n = 2 through rid and (inner 
loop) m = 1 through n — 1 

4.4 

/ s 

(sec </>Pn,m) = [(2 n- 1) sin cj) A„_ liTO (n, m)(sec<^P„_ lim ) 

- (n + m- 1) A n _ 2 , m (n,m)(sec</>P„_ 2 ,m)] — - — (4.9) 

v v ' n — m 

4.5 

where (sec</> P„_i,„) = 0 for n = 1 through n^. 

Recursion for Sectorial ALF Derivatives (cos 4>P' n ,n) 

This recursion has no normalization parameters because the input and output 
are the same degree and order. For n = 1 through rid 

(cost j>P' n , n ) = -rising (sec (j>P n , n ) (4-10) 
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Recursion for Tesseral ALF Derivatives (cos (t>P' n ,m) 

This recursion utilizes Parameter 4.4. For n = 2 through rid and (inner loop) 
to = 1 through n — 1 

(cos (jrP'n^m) = — n sin </> (sec </> P„,m) 

+ (n + m) A„_i im (?r,TO)(sec^P„_i iTra ) (4.11) 

4.4 

4.2.3 Example of Normalized Lear Recursions 

This section analytically demonstrates using normalized recursion relationships 
to generate Legendre polynomials and ALFs for a gravity model with a degree 
and order of four. Each of the final answers are written first simplified and 
then with its normalization factored out to demonstrate equivalence with its 
unnormalized value. These can be found in the examples outlined in Ref. [2]. 


Zonal Legendre Polynomials 


Pi 


h 


P4 


- [3sin(/> A„_i(2)Pi — A„_ 2 (2)P 0 ] 


3 sin (j) y ^ a/ 3 sin (j) — a/5 
^ a/ 5 sin 2 </> — ^ a/5 = a/ 5 sin 2 </> — 

[5sin^.A„_ 1 (3)P 2 -2A„_ 2 (3)P 1 ] 

5 sin (f> y ^ ^ - a/ 5 sin 2 <j> — - a/5^ — 2 a/ 3 sin (f> 


5-3 5 

— — a/7 sin 3 (f> — - a/7 sin cj) — 2\fl sin </> 


/ a/7 sin 3 (j) — ^ a/7 sin (f> = a/7 ^ / sin 3 4> — 3 sin <fi 


- [7sin</»A„_ 1 (4)P 3 -3A„_ 2 (4)P 2 


7 sin^y| ^\/7 sin 3 $ — ^\/7 sin<^ 

- 3 /I 
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1 [7-5-3 

4 [ 2 

105 • 4 , 

= sm 0 

O ” 


sin 4 (j) 


7-3-3. 2 , 3-3-3. 2 , 3-3 

— ^ — sm 0 — sm 0 + — 




3 

8 


Zonal Legendre Polynomial Derivatives 


P '2 = A n _ 1 (2)[sin0P' 1 +2P 1 ] 
5 

3 L~ 

= 3A/5sin0 = a/ 5 (3sin0) 


sin 0 a/ 3 + 2a/ 3 sin 0 


P's = A„_i(3) [sin0P' 2 + 3P 2 ] 


sin 0 ^3v/5 sin cfij + 3 a/ 5 sin 2 0 — ^ v/5 


7 

5 _ 

3v/7 sin 2 0 + vT sin 2 0 — ^ a/7 


^ a/ 7 sin 2 0 — / a/7 = a/7 f ^ sin” 0 — - 


15 . 2 

— sin z 
2 


P’a = 


A„-i(4) [sin 0P' 3 + 4P 3 ] 

\J ^ sin0 a/ 7 sin 2 0 — ^ a/7^ + 4 ^a/ 7 sin 3 0 — ^\/7 sir 


3-15 


3-3 


3-4-5 


4-3-3 


sm 0 sm ( 

2 Y 2 


■ sm 0 — 


smi 


105 . o 45 . , / 35 , , 15 . 

sin 0 sm 0 = 3 — sin 0 sm < 

2 2 \2 2 


Sectorial (Diagonal) ALFs 


(sec 0P 2 , 2 ) 


(sec 0P 3 , 3 ) 


3 cos 0 A„_i >n _i (2) (sec 0 P M ) 

3cos 0^y| v/ 3 

1 l/5^ 
-a/15cos0 = - \ - (3 cos 0) 

2 2 V 3 v ' 

5 cos 0 A„_i >n _i (3) (sec 0 P 2j2 ) 
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= 5cos<^y|^\/T5cos<^ 

= ^/f COsV= ^ (l5COsV) 


(sec 0 P 4i 4 ) = 7cos0Ara-i,n-i(4)(sec<£P 3i3 ) 

I/ 9 / 1/35 2 \ 

= 7co ^7V8 UVT COS 0 


= ^V35cos 3 8>= 1^4(105cos 3 8>) 


Tesseral ALFs: Row 2 (n = 2) 


(sec0P 2a ) = 3sin^A n _i !m (2, l)(sec^Pi i i) -2A„_ 2>m (2, l)(sec^P 0 ,i) 
= 3sm(f>\ - — — >/3 — 2 y/O ■ 0 

V 0*0 


= Vl5 sin (j> = (3 sin </>) 


Tesseral ALFs: Row 3 {n = 3) 


(sec^P 3 ,i) = [5sin<^A n _i >m (3, l)(sec0P 2 ,i) - 3A„_ 2 , m (3, l)(sec^P l> i)] - 

= 2 5sin ^ \/M (^ sin ^) ^ 3 \/ 4 ^ 3 ^ 

5 Fu3 . „ , 1 /7^3 


5 7 • 3 . 2 , 1 /7 - 3 

JV / — S m 

5 /21 . 2 , 1/2! /7/15.2, 3 

2 V 2 2V2 V 6 V 2 2 


(sec ^P 3 , 2 ) = 5sin<£A„_i jm (3, 2)(sec0P 2 , 2 ) - 4A n _ 2 , m (3, 2)(sec</>Pi i2 ) 


= 5 sin 


0 ^vT5cos</)^ -4\/0-0 


5 / 7 • 5 • 3 
2 V 5-5 


sin <f> cos 


^ 2 / 7 

= — v 7 105 sin (j> cos (f> = -\/— (15 sin cos (^>) 

Z Z y -Lo 
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Tesseral ALFs: Row 4 (n = 4) 

The example in the Lear document which corresponds to the following example 
contains a typographical error. The final result for the ALF should be 


, , 35 . 3 15 . 

(sec 0 P 4; i) = — sin 0 — — sin < 


(sec0P 4>1 ) = [7sin0A n „i im (4, l)(sec0P 3j i) - 4A„_ 2>m (4, l)(sec0P 2 ,i). 


7 . 

- sm< 
3 


3-9/5 / 21 


1 / 21 


5 • 7 1 2 V 2 Sm 0 2 V 2 


4 /3 - 2 - 9 


3 V 5-4-5 


15 sin 


7-5 /3 -9-7-3 


7 /3 • 9 • 7 • 3 


3-2 V 5-7-2 


■ sm 0 — 


3- 2 V 5-7-2 


sin 0 


4 /3 • 2 • 9 • 3 • 5 


5-4-5 


sin 0 


21 /5 . 3 , 9/5 


TV2 Sm 2 V 2 Sm ' 

„ /T /35 . 3 , 15 . 

= 3 VioU sin 


(sec 0P 4 , 2 ) = [7sin0A„_i irn (4,2)(sec0P 2i 2) - 5A„_ 2>m (4, 2)(sec0P 2j2 ) 


7 / 2*9/1 \ 

2 sin 0 Y gTy ( 2V / iO5sin0cos0J 

2 • 9 /I 

-v 15 cos 0 


2 V 6-5-5 V2 
7 /2 - 3 - 3 - 3 - 5 - 7 


5 /2 • 3 • 5 • 3 • 3 


2-3-7 


- sm m cos 0 

4 


2 • 3 • 5 • 5 


COS 0 


— a/5 sin 2 0 cos 0 — - a/5 cos 0 
4 4 


1 /I / 105 


2 V 5 V 2 ™-^ c °s0-L c os, 


(sec 0P 4 , 3 ) = 7sin0A„_i, m (4,3)(sec0Pi,i) - 6A„_ 2jm (4, 3)(sec0P 2 , 3 ) 


= 7sin0\/y^ | iyycos 2 0j -6/0-0 


7 /3 -3-5-7 


2-7-7 


• sm 0 cos 
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= s\[r sin ^ cos2 ^ = ^\/^( 105sin ^ cos2 ^) 

Sectorial (Diagonal) ALF Derivatives 


(cos0P'i,i) = 


— sin0(sec^)Pi ) i) 

— a/3 sm4> = v / 3(— sin^) 


(cos 4> P' 2 , 2 ) = —2 sin <j> (sec <f> P 2i2 ) 

= — 2sin(/> ^a/15cos(^ 

1 [b 

= -Vl 5 sm(pcos(j>=-\-(— 6 sin(f>cos(j)) 

Zi \ o 


(COS (f) P' 3 , 3 ) = 


— 3sin</> (sec^P 3i3 ) 

, • , f 1 /35 2 ^ 

—3 sm 0 I 9 \^ 2 COS V 


3 / 35 


sin 0 cos 


(cos (j) P' 4 ,, 4 ) = 


2 V 2 
— 4 sin </> (sec <)) p 4,4 


4 > = ^ y (-45 sin (f> cos 2 </>) 


= —4 sin (j) ^\/35cos 3 

3 1 / 1 

= — - a/35 sin (f> cos 3 </> = g y 35 (“420 sin (j) cos 3 </>) 

Tesseral ALF Derivatives: Row 2 (n = 2) 


(cos^P' 2 ,i) = -2sin^(sec0P 2 ,i) + 3A„_i >m (2, l)(sec^Pi,i) 


= —2 sin cj) ( \/T5 sin /’) + 3 \/ o o v/ ^ 

o ' o 


= — 2A/l5sin 2 (j) + a/ 15 = y/jj (— 6 sin 2 (j) + 3) 

Tesseral ALF Derivatives: Row 3 (n = 3) 


(cos</>P' 3 ,i) 


—3 sin (/> (sec </> P 3) i ) + 4A n _i, m (3, l)(sec</>P 2 ,i) 
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= — 3sin< 


wh 


15 sin 


15 / 21 . , 3 21 . , r- „ , - . 

- y Y y sin 3 (j) + - J — sin <t> + V2 • 7 • 4 • 3 sin < 


15 21 


11 21 


sm 


45 


2 V 2 

, 33 


smi 


= \ — \ sin 

6 V 2 


■ sm ( 


(cos(/> P'3,2) = -3sin^(sec</>P 3j2 ) + 5A„_i, m (3, 2)(sec </>P 2 , 2 ) 

= — 3sin</> ^ \/l05 sin (f> cos (j^J + 5y^ ^ \/l5 cos (j> 


3 1 

— -%/lC)5sin 2 (f>cos(J)+ -\/l05cos<^ 

1 

= - \ — (—45 sin 2 q i cos <j> + 15 cos (/>) 
Z V -L 0 


Tesseral ALF Derivatives: Row 4 (?r = 4) 

The example in the Lear document which corresponds to the following example 
contains a typographical error. The final result for the ALF should be 

135 15 

(cos (j> P' 4 , 1 ) = —70 sin 4 <j> H — — sin 2 0 — — 


(cos^P'4,1) = -4sin^(sec0P 4 ,i) + 5A„_i, m (4, l)(sec0P 3i i) 

. . , ^21 [5 . 3 , 9 [5 . \ 


_ 73-9/5 /21 . 2 , 1/21 

Sl, w ! VT™ ^VT 


— 21\/l0sin 4 + 9VlO sin 2 0 + ^4 i/L sin 2 ^ 


5 /3 • 9 • 7 • 3 
2 V 2-5-7 


- 21 \/l 0 sin 4 0 + y ^ sin 2 0 - | 


= 3\/ y ( —70 sin 4 - 


135 


15 


sm 0 

2 ^2 


(cos (j) P' 4j2 ) = -4sin^)(sec^P 4i2 ) + 6A n _!, m (4, 2)(sec^P 3 , 2 ) 
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( 21 3 \ 

— a/5 sin 2 (f> cos <f> — - a/5 cos </> j 

+ 6\l - — - ( -a/105 sin 0 cos 0 
V 6 • 7 \2 y 


= — 21\/5 sin 3 4> cos <^> + 3 a/ 5 sinc/i cos </> 

„ / 9- 3- 5- 7 . , 

+ 3 W sin d> cos </) 

V 3-7 

= — 21-\/5sin 3 tficoscj) + 12y/5 sin <j> cos <j) 

= ^ (—210 sin 3 (j> cos (j) + 120 sin (f> cos (/>) 

(cos (f) PA, 3 ) = -4 sin ^ (sec <^ 4 , 3 ) + 7A„_i >m (4, 3)(sec </> P 3j3 ) 


= — 4sin< 


3 / 35 


sm m cos 


+ 7 


1 / 35 


7 • 7 V 2 


cos 


— 3a/ 70 sin 2 <f> cos 2 </> + ^ \j^~ cos2 & 

= ^ (—420 sin 2 <f> cos 2 </> + 105 cos 2 </>) 


4.3 Gottlieb Algorithm 

4.3.1 Basis of the Gottlieb Approach 

Mueller [10] developed an efficient algorithm for solving the gravitational po- 
tential function which, like Pines, defined special functions by reallocating the 
factors in each term of the potential function. Gottlieb [3] then defined the 
gradient of Mueller’s potential function in terms of the partial derivatives with 
respect to these special functions. Recursions for the partial derivatives of the 
special functions were then developed to resolve the gravitational acceleration. 

4.3.2 Normalized Gottlieb Algorithm 

Each of the recursion equations of the Gottlieb algorithm is now normalized 
using the normalization parameters of Section 3.4. Each of the original equa- 
tions can be found in Ref. [3]. Only the normalized equations are presented 
here and should replace their analogues in the original algorithm. The added 
normalization parameters are denoted with an underbrace or overbrace to bring 
attention to the changes from the equations in the original document. 

Normalizing the Gottlieb recursions requires a few new normalization pa- 
rameters in addition to the parameters derived for a normalized Lear imple- 
mentation. The new normalization parameters are outlined in Table 4.2. In 
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^n,m+l TYl) — 

Hn,m 

N n ,m-\- 1 


^n,m+l(^5 0) — 

ii 

O t-H 

!-l,m-l(n,ra) = 

Nn,m 

H n—l,m 



(n + to + l)(n — m)( 2 — d 0 ,m) 


(2n + 1)(2 — So, m ) 


(2 n — l)(n + m)(n + to — 1)(2 — <5o,m-i) 


(4.12) 

(4.13) 

(4.14) 

(4.15) 


Table 4.2: Normalization parameters A needed for a normalized Gottlieb algo- 
rithm 


this section, sin </> has been replaced by e to keep a consistent notation with the 
original document. 

Recursion for Zonal Legendre Polynomials P° = P n 

This recursion utilizes Parameters 4.1 and 4.2. For n — 2 through rid 

= P n = 1[(2 n - l)e A„_i(n) P„-i - (n - 1) A„_ 2 (n) P„_ 2 ] (4.16) 

4.1 4.2 

where Pq = 1 and Pi = NiP\ = v^e- 


Recursion for Sectorial and Tesseral ALFs P™ 


This recursion utilizes Parameters 4.5 and 4.14. For n = 2 through rid and 
(inner loop) m = 1 through rid 


P™ = A n _ 2 >m (n,TO)P ™_ 2 + (2n 

s v y 

4.5 


1) A„_ ljm _i (n, to) P™_Ti 

\ ^ / 

4.14 


(4.17) 


where P/ = 1 . 


Recursion for Intermediate Sum H n 

This recursion utilizes Parameters 4.13 and 4.12. For n = 2 through rid 


4.13 


Hn — C n , 0 ^n,m+l(^5 0) H n 


n pm + 1 

^n,m+l (^5 ~m. H” *$n,ra*$m) (4.18) 


m=l 


r" 


4.12 
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Chapter 5 

Verifications and 
Conclusions 


Each of the algorithms addressed in this paper were at first coded directly 
from the equations or code provided in the original sources and then tested 
against normalized implementations of each. Surprisingly, the results showed a 
fairly large discrepancy between the results of the normalized and unnormalized 
implementations for two of the three algorithms. Through numerical analysis 
and extremely scrutinized debugging, it became evident that the generator for 
the ALFs was the source of error. 

Since the normalized Pines implementation was provided to the authors al- 
ready coded but its results appeared more stable than the unnormalized code, 
it was suspected that the two codes were in some way different. As an exper- 
iment, the unnormalized Pines code was replaced by a copy of the normalized 
version of the code which had been “unnornralized,” revealing that the ALF 
generator was in fact different between the two Pines implementations utilized 
in the first test. Further experimentation verified that the stability of the al- 
gorithms depended largely on the stability of the ALF generator used. Finally, 
the results of all testing appeared to indicate that normalization amplifies any 
inherent noise and error in each of the algorithms, a conclusion which further 
drove the development of additional conclusions and recommendations. 


5.1 Preliminary Test Conditions 

To test both the normalized and unnormalized implementations of each algo- 
rithm (six total) for agreement, acceleration vectors were computed at a set of 
positions around a given central body. The Moon was chosen for the central 
body, and the LP150Q mass coefficients were utilized for the test. The mass 
coefficients were unnornralized using the recursions in Section 3.2 to pass to 
the unnormalized algorithms. Tests utilized central-body- fixed position vectors 
with latitudes ranging from —90° to +90° in 30° increments and at each longi- 

25 
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tude from —150° to +180° in 30° increments. Each position vector was given a 
200-kilometer altitude above the lunar reference radius of 1738 kilometers. 

At each of these positions, the acceleration was computed with all six al- 
gorithms with a variety of gravity model sizes. As a first check to ensure the 
algorithms functioned properly, a 0 x 0 model was tested to obtain the central 
body acceleration. Square models 2x2 through 50 x 50 were then tested, fol- 
lowed by the non-square models 50 x 0 through 50 x 49. Finally, the “extreme” 
cases of 125 x 125 and 150 x 150 were tested to ensure the normalized models 
in fact converged at relatively high degrees and orders. 

The DeMars implementation of normalized Pines was considered the baseline 
model because it was based on the extensive stability studies of Lundberg and 
Schutz [5] . Error, defined for each tested iteration as the magnitude of the delta 
vector between the DeMars-calculated acceleration vector and the calculated 
accleration vectors from each of the other algorithms, was considered acceptable 
if the order of magnitude was 10~ 18 or smaller. This is the order of magnitude 
of ten times machine epsilon for the acceleration magnitudes tested, which was 
obtained by passing various acceleration vector magnitudes from the DeMars 
subroutine as arguments to the MATLAB eps function. 


5.2 Preliminary Results 

The output of the first MATLAB script to test the algorithms is listed in Ap- 
pendix A. The preliminary results showed that Lear (both normalized and 
unnormalized) and Pines unnormalized had negligible error for all gravity mod- 
els and locations around the central body. Gottlieb had small but noticeable 
error in equatorial locations in the square models 40 x 40 and larger. All of 
the Gottlieb non-square models had large but consistent error everywhere. The 
error became especially pronounced in the large 125 x 125 and 150 x 150 models 
at the equator. These findings necessitated further analysis of the error which 
developed at the equator. 


5.3 Further Analysis 

The tests outlined in Section 5.1 were modified to only calculate accelerations 
with all the square models from 2x2 through 150 x 150. This would allow a 
trend to emerge when plotting the error versus the degree and order of the model, 
as shown in Figure 5.1. At the equator, the error in both the normalized and 
unnormalized Gottlieb models grows slowly with increasing degree and order and 
suddenly diverges to positive infinity when the degree and order approaches 150. 
The rate of growth of error in Gottlieb is drastically larger as seen by comparing 
the scales of the y-axis of Figures 5.1 and 5.2. The beginning of the divergence 
of Gottlieb from the other two models, which remain closely in line with each 
other, in the unstable implementations can be clearly seen in Figure 5.3. 

As an experiment, this same test was run using a known unstable ALF 
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Figure 5.1: Acceleration error magnitude for unnormalized models at equator 
(0 = 0 °) 
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Figure 5.2: Acceleration error magnitude for normalized models at equator 
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Figure 5.3: Acceleration error magnitude for unnormalized models at equator 
(</> = 0°) with degree and order 60-85 


generator 1 in the Pines algorithm. The error behavior of this known unstable 
algorithm at the poles mirrored the behavior of the Gottlieb error at the equator, 
as shown in Figure 5.4. 

Normalized implementations of the unstable Pines showed that the error 
was much larger than the error of the unnormalized unstable Pines, as seen in 
Figure 5.5. This large disparity between normalized and unnormalized error 
resembled the error of the two Gottlieb implementations as well. 

5.4 Conclusions 

Four primary conclusions can be drawn from the data presented in this paper: 

• Pines (as implemented by DeMars) and Lear algorithms are sta- 
ble because they use a stable ALF recursion. It is worth noting 
that virtually the same recursion equation [4, Recursion I] is used for gen- 
erating ALFs in the unnornralized implementations of the Lear and Pines 
(DeMars) algorithms. The very similar behavior between the two algo- 
rithms can thus be explained by the similarity in their ALF generators. 

• Gottlieb and Pines algorithms, as originally published, are un- 
stable due to unstable ALF recursions. The apparently unstable be- 
havior of the Gottlieb algorithm is presumed to be the result of an unstable 
ALF generator used in the algorithm. This conclusion is motivated by the 
similar signature of the error data with a known unstable ALF generator, 
the Pines algorithm with Unnornralized Recursion IV from Lundberg and 

1 Unnormalized Recursion IV from Lundberg and Schutz [4], 
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Figure 5.4: Acceleration error magnitude for unnormalized models (unstable 
Pines) at south pole (0 = —90°) 
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Figure 5.5: Acceleration error magnitude for normalized models (unstable 
Pines) at south pole (</> = —90°) 
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Schutz [4] with which it was originally published. The error is location 
specific, much like the error that develops in the inherently unstable Pines 
implementation, albeit in a different location. The fact that both error 
signatures are latitude-dependent implies the ALF generator is the source 
of the instability, since the argument of ALFs used for spherical harmonic 
expansion is always sm<t>. 

• Normalization of recursions amplifies numerical noise. Assuming 
Gottlieb is, in fact, using an unstable ALF generator, it would appear that 
normalization of unstable ALF generators increases this instability and 
amplifies the error dramatically. This conclusion is supported by observing 
the same amplified error in both a known unstable ALF generator in Pines 
and the presumed-unstable ALF generator in Gottlieb in their normalized 
implementations relative to their unnormalized equivalents. 

• Unnormalized algorithms provide perfectly valid results at high 
degree and order as long as coefficients can be reliably unnor- 
malized. Normalized and unnormalized implementations of Pines as well 
as Lear algorithms agree very well with each other. This leads to the con- 
clusion that it is safe to use unnormalized algorithms as long as proper 
unnormalization of the coefficients is performed, such as using the recur- 
sions in Section 3.2. If a relatively small gravity model is always desired, 
such as in current Mission Control Center software where the Earth gravity 
model is traditionally limited to degree and order 7, it is perfectly accept- 
able to continue the practice of implementing unnormalized algorithms for 
calculating gravitational acceleration. 

5.5 Recommendations 

Two recommendations are presented by the authors: 

• Gottlieb and Pines algorithms should not be implemented di- 
rectly as published. Unless a more stable ALF generation scheme is 
implemented, it is recommended that the Gottlieb algorithm be omitted 
from any implementation which incorporates large models due to the po- 
tential for instability. Developing an improved Gottlieb algorithm which 
implements a stable ALF generator is left for future work. 

The Pines algorithm has been stabilized by Lundberg and Schutz [4], so 
their recursions should be implemented with a Pines algorithm instead of 
the ALF generator used in the original Pines paper. 

• Normalized implementations are better suited to software pack- 
ages than unnormalized algorithms. Normalized algorithms are not 
only better suited for acceleration calculation with larger gravity mod- 
els, their consistency with unnormalized algorithms makes them desirable 
for implementation in software packages seeking to maintain versatility 
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and robustness when computing gravitational acceleration. Normalized 
algorithms, when properly implemented, should always return a valid ac- 
celeration given any set of normalized coefficients and a valid degree and 
order. 


This document has been reviewed for Proprietary, SBU, and Export Control (ITAR/EAR) and has been determined to be nonsensitive. 
It has been released to the public via the NASA Scientific and Technical Information (STI) Process DAA 23097. 



32 CHAPTER 5. VERIFICATIONS AND CONCLUSIONS 


THIS PAGE IS INTENTIONALLY LEFT BLANK. 


This document has been reviewed for Proprietary, SBU, and Export Control (ITAR/EAR) and has been determined to be nonsensitive. 
It has been released to the public via the NASA Scientific and Technical Information (STI) Process DAA 23097. 



Appendix A 

Preliminary Results 


This section contains the output of the initial MATLAB test script. These ini- 
tial tests of the algorithms were used to identify test cases which needed further 
analysis. Cases with unusually large error are boxed like this . Note: Error is 


defined as the magnitude of the delta vector between each algorithm’s acceler- 
ation vector and the Normalized Pines acceleration vector. See Section 5.1 for 
additional information. 


Maximum central body unnormalized Pines error: 2.65574e-019 
at lat/lon: -30/-60 
degree x order: 0x0 

Maximum central body unnormalized Lear error: 4.8487e-019 
at lat/lon: 0/-120 
degree x order: 0x0 

Maximum central body unnormalized Gottlieb error: 6.50521e-019 
at lat/lon: -90/-150 
degree x order: 0x0 

Maximum central body normalized Lear error: 0 
at lat/lon: -90/-150 
degree x order: 0x0 

Maximum central body normalized Gottlieb error: 0 
at lat/lon: -90/-150 
degree x order: 0x0 


Maximum square unnormalized Lear error: 1.99033e-018 
at lat/lon: 0/30 
degree x order: 48x48 


Maximum square unnormalized Gottlieb error: 2.7959e-018 
at lat/lon: 30/-150 
degree x order: 49x49 


Maximum non-square unnormalized Lear error: 2.28713e-018 
at lat/lon: 30/-90 
degree x order: 50x9 
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APPENDIX A. PRELIMINARY RESULTS 


Maximum non-square unnormalized Gottlieb error: 8.29053e-008 


at lat/lon: 60/-150 
degree x order: 50x0 

Maximum square Pines error: 2.65574e-019 
at lat/lon: -30/-150 
degree x order: lxl 

Maximum square Lear error: 2.48422e-019 
at lat/lon: 60/60 
degree x order: 46x46 


Maximum square Gottlieb error: 1.32303e-016 


at lat/lon: 0/-150 
degree x order: 49x49 

Maximum non-square Pines error: 4.96844e-019 
at lat/lon: -60/-60 
degree x order: 50x13 

Maximum non-square Lear error: 2.65574e-019 
at lat/lon: -30/150 
degree x order: 50x25 


Maximum non-square Gottlieb error: 1.43628e-016 


at lat/lon: 0/-150 
degree x order: 50x26 

Maximum square normalized Lear error 125: 9.00606e-019 
at lat/lon: 30/-90 
degree x order: 125x125 


Maximum square normalized Gottlieb error 125: 5.96756e-005 


at lat/lon: 0/-120 
degree x order: 125x125 

Maximum square normalized Lear error 150: 2.7959e-018 
at lat/lon: 30/150 
degree x order: 150x150 


Maximum square normalized Gottlieb error 150: 1.17375 


at lat/lon: 0/-30 
degree x order: 150x150 
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